#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>

#define ROOT 0

int main(int argc, char *argv[]) {
    int J_ORDER, J_ROW_TEST, J_ITE_MAX, i= 0, j= 0, k, l, iteracoes= 1, criterio_de_convergencia = 0, total, rank, size, parte, flag = 1;
    float J_ERROR, erro= 0, maior_X, aux= 0, aux2;

    FILE *f = fopen(argv[1], "r");
    if (f == NULL)
        return 0;

    fscanf(f, "%d %d %f %d", &J_ORDER, &J_ROW_TEST, &J_ERROR, &J_ITE_MAX);

    float *A, *A_, *auxMatriz, B[J_ORDER], B_[J_ORDER], X[2][J_ORDER], multiplic[J_ORDER], resultado[J_ORDER], somatoria[J_ORDER];

    A = (float *)malloc(J_ORDER* J_ORDER* sizeof(float));
    A_ = (float *)malloc(J_ORDER* J_ORDER* sizeof(float));
    auxMatriz = (float *)malloc(J_ORDER* J_ORDER* sizeof(float));

    double inicio = MPI_Wtime();

    MPI_Init( &argc,&argv);
    MPI_Comm_rank( MPI_COMM_WORLD, &rank);
    MPI_Comm_size( MPI_COMM_WORLD, &size);

    parte = J_ORDER / size;

    if (rank == 0) {
        /* leitura e preparação das matrizes */
        for((i= 0); (i< J_ORDER); (i++)) {
            for((j= (i* J_ORDER)); (j< (J_ORDER* (i+ 1))); (j++)) {
                fscanf(f, "%f", &A[j]);
                A_[j] = A[j];
            }
        }
        for (i = 0; i < J_ORDER; i++) {
            fscanf(f, "%f", &B[i]);
            B_[i] = B[i];
        }

        j = 0;
        #pragma omp parallel for
        for (i = 0; i < J_ORDER; i++) {
            j = i*(J_ORDER+ 1);
            B_[i] /= A[j];
            X[0][i] = B_[i];
        }
        #pragma omp parallel for shared(A_, A) private (i, j)
        for(i= 0; i< J_ORDER; i++) {
            for(j= (i* J_ORDER); j< (J_ORDER* (i+ 1)); j++) {
                if (j!= i*(J_ORDER+ 1)) A_[j] /= A[i*(J_ORDER+ 1)];
            }
            A_[i*(J_ORDER+ 1)] = 0;
        }

        MPI_Bcast(X[0], J_ORDER, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
        /* envias as partes da A_ */
        MPI_Scatter(&(A_[0]), parte * J_ORDER, MPI_FLOAT, &(auxMatriz[0]), parte * J_ORDER, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
        for (k = 0; k < J_ITE_MAX && !criterio_de_convergencia; k++) {

            for((i= 0); (i< parte); (i++)) {
                resultado[i] = 0;

                for((j= (i* J_ORDER), l= 0); (j< (J_ORDER* (i+ 1))); (j++, l++)) {
                    resultado[i] += auxMatriz[j] * X[0][l];
                }
            }
            /* reúne os calculosdas somatorias */
            MPI_Gather(resultado, parte, MPI_FLOAT, somatoria, parte, MPI_FLOAT, ROOT, MPI_COMM_WORLD);

            /* X[k+1] */
            for (i = 0; i < J_ORDER; i++)
                X[1][i] = B_[i] - somatoria[i];

            aux = 0;
            for (i = J_ROW_TEST* J_ORDER, j= 0; i < (J_ROW_TEST+ 1)* J_ORDER; i++, j++)
                aux += A[i] * X[1][j];

            for (i = 0; i < J_ORDER; i++)
                X[0][i] = X[1][i];

            iteracoes++;

            MPI_Bcast(X[0], J_ORDER, MPI_FLOAT, ROOT, MPI_COMM_WORLD);

            if (fabs(B[J_ROW_TEST] - aux) < J_ERROR)
                criterio_de_convergencia = 1;
            MPI_Bcast(&criterio_de_convergencia, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
        }
        printf("\n---------------------------------------------------------\nIterations: %d\nRow Test: %d => [%f] =? %f\n---------------------------------------------------------\nTempo: %lf segundos\n\n", iteracoes-1, J_ROW_TEST, aux, B[J_ROW_TEST], (MPI_Wtime() - inicio));
    }
    else {
        MPI_Bcast(X[0], J_ORDER, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
        MPI_Scatter(&(A_[0]), parte * J_ORDER, MPI_FLOAT, &(auxMatriz[0]), parte * J_ORDER, MPI_FLOAT, ROOT, MPI_COMM_WORLD);

        for (k = 0; k < J_ITE_MAX && !criterio_de_convergencia; k++) {
            /* A_ x B_ */
            for((i= 0); (i< parte); (i++)) {
                resultado[i] = 0;

                for((j= (i* J_ORDER), l= 0); (j< (J_ORDER* (i+ 1))); (j++, l++)) {
                    resultado[i] += auxMatriz[j] * X[0][l];
                }
            }
            MPI_Gather(resultado, parte, MPI_FLOAT, somatoria, parte, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
            MPI_Bcast(X[0], J_ORDER, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
            MPI_Bcast(&criterio_de_convergencia, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
        }
    }
    MPI_Finalize();
    return 0;
}
